#include "nnec/nnec_utils.h"
#include "nnec/nnec_errors.h"

using namespace std;

namespace NNEC
{
    int integration_romberg(INTEGRAND *integrand, double a, double b, int n_max,
                            double epsabs, double epsrel, double *result, int *eval_count)
    {
        *result = 0.0;

        if(n_max < 1 || n_max > 30) {
            return NNEC_ERROR_INVALID;
        }
        if(epsabs < 0.0 || epsrel < 0.0) {
            return NNEC_ERROR_BADTOLERANCE;
        }

        auto previous_result{make_unique<double[]>(n_max)};
        auto current_result{make_unique<double[]>(n_max)};
        if(!previous_result || !current_result) {
            return NNEC_ERROR_ALLOCATION;
        }

        double h = 0.5 * (b - a);

        previous_result[0] = h * (integrand->function(a, integrand->parameter) + integrand->function(b, integrand->parameter));
        *eval_count = 2;

        for(int i = 1; i < n_max; i++) {
            double sum = 0.0;
            int two_i = 1 << i;

            for(int j = 1; j < two_i; j = j + 2) {
                sum = sum + integrand->function(a + j * h, integrand->parameter);
                (*eval_count)++;
            }

            current_result[0] = sum * h + 0.5 * previous_result[0];

            double four_j = 4.0;
            for (int j = 1; j <= i; j++)
            {
                current_result[j] = (four_j * current_result[j - 1] - previous_result[j - 1]) / (four_j - 1.0);
                four_j = 4.0 * four_j;
            }

            double err = abs(current_result[i] - previous_result[i - 1]);
            if((err < epsabs) || (err < epsrel * abs(current_result[i])))
            {
                *result = current_result[i];
                return NNEC_SUCCESS;
            }

            swap(previous_result, current_result);

            h = 0.5 * h;
        }

        *result = previous_result[n_max - 1];

        return NNEC_ERROR_MAXITERATION;
    }
}
